### LCV source file


commarobust <- function(fit){
  require(lmtest)
  require(sandwich)
  coeftest(fit,vcovHC(fit,type = "HC2"))
}

getrobustses <- function(fit){
  robustfit <- commarobust(fit)
  #N <- length(fit$residuals)
  #p <- 1/N
  returnvec <- robustfit[,2]
  returnvec[is.nan(returnvec)] <- .00004
  #returnvec <- round(returnvec, 4)
  return(returnvec)
}


getcoefs <- function(fit){
  robustfit <- commarobust(fit)
  returnvec <- robustfit[,1]
  returnvec <- round(returnvec, 6)
  return(returnvec)
}

makerobustseslist <- function(fitlist){
  return(lapply(fitlist, FUN=getrobustses) )
}
makecoefslist <- function(fitlist){
  return(lapply(fitlist, FUN=getcoefs) )
}



# function approximating true RA protocol
strata_ra <- function(stratum, treat){
  treat <- as.character(treat)
  assign <- rep(NA, length(stratum))
  strata <- unique(stratum)
  for(i in 1:length(strata)){
    n_stratum <- sum(stratum==strata[i])
    n_pub <- sum(treat[stratum==strata[i]]=="pub")
    n_fol <- sum(treat[stratum==strata[i]]=="fol")
    n_org <- sum(treat[stratum==strata[i]]=="org")
    rand <- sample(1:n_stratum, replace=FALSE)
    assign_stratum <- rep(NA, n_stratum)
    assign_stratum[rand[1:n_pub]] <- "pub"
    assign_stratum[rand[(1+n_pub):(n_pub+n_fol)]] <- "fol"
    assign_stratum[rand[(1+n_pub+n_fol):n_stratum]] <- "ord"
    assign[stratum==strata[i]] <- assign_stratum
  }
  return(assign)
}


# helper functions
se_mean <- function(x){
  x.nona <- x[!is.na(x)]
  n <- length(x.nona)
  return(sd(x.nona)/(sqrt(n)))
}

f_tester <- function(x,z){
  fit <- lm(x ~ z)
  if(sum(fit$residuals^2) < 1e-10) {return(p=1)} # dealing with perfect fits
  f_stat <- summary(fit)$fstatistic
  return(p=pf(q = f_stat[1], df1 = f_stat[2],df2 = f_stat[3],lower.tail = FALSE))
}

add_parens <- function(x, digits=3){
  x <- as.numeric(x)
  return(paste0("(", sprintf(paste0("%.", digits, "f"), x), ")"))
}

format_num <- function(x, digits=3){
  x <- as.numeric(x)
  return(paste0(sprintf(paste0("%.", digits, "f"), x)))
}



